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ABSTRACT 

Brief  descriptions  are  presented  of  several  of  IIT  Research 
Institute's  analyses  and  associated  computer  codes  relating  to 
shock  wave  propagation  through  earth  media  and  to  the  response 
of  buried  structures  to  ground  shock  loading. 
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I.  INTRODUCTION 


This  report  contains  a compilation  of  summaries  of  anal- 
yses developed  at  IIT  Research  Institute  pertaining  to  wave 
propagation  through  earth  media  and  soil-structure  interaction 
as  a result  of  nuclear  detonation  or  other  high  intensity  wave 
sources . 

The  summaries  can  basically  be  divided  into  two  groups 
The  first  group  describes  analyses  which  were  concerned  with 
determination  of  the  free-field  behavior  of  earth  media  sub- 
jected to  shock  wave  propagation  through  it.  The  second  group 
is  composed  of  studies  devoted  to  the  prediction  of  horizontal 
and/or  vertical  response  of  buried  tunnel  and  missile-silo 
structures  to  ground  shock  loading. 

References  for  each  analysis  are  cited  by  number  immedi- 
ately following  the  title.  The  corresponding  reference  infor- 
mation is  contained  in  the  list  of  references  at  the  end  of 
the  report. 

II.  FREE-FIELD  ENVIRONMENT 

The  studies  associated  with  determination  of  the  free- 
field  behavior  of  the  earth  media  have  been  separated  either  into 
a category  in  which  the  analyses  are  based  on  a hydrodynamic 
description  of  material  behavior,  or  one  in  which  the  analyses 
characterize  the  soil  as  a solid.  There  is  no  essential  differ- 
ence between  the  two  approaches  for  one-dimensional  problems. 

A . Hydrodynamic  Models 

1 . SW1MM  Code  (Ref.  1) 

The  SWIMM  (Stress  Waves  in  Multilayered  Media)  Code  was 
written  in  FORTRAN  II  for  use  on  an  IBM  7094  digital  computer. 

The  purpose  of  the  code  is  to  calculate  one-dimensional  uniaxial 
stress  wave  attenuation  in  layered  media.  The  main  uses  thus 
far  of  the  code  have  been  in  accomplishing  calculations  involved 
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in  the  analysis  of  various  blast  shield  designs  for  use  with 
a nuclear  reactor  in  the  advent  of  a malfunction  due  to  over- 
heating. However,  the  code  can  also  be  directly  applied  to 
predict  the  air-induced  ground  shock  environment  created  by 
an  H.E.  or  nuclear  explosion. 

The  media  may  contain  from  one  to  forty  separate  layers, 
as  illustrated  in  Fig.  1,  with  a prescribed  loading  at  the  top 
surface  and  a base  boundary  condition  corresponding  either  to 
a rigid  wall  (i.e.,  velocity  equal  to  zero),  or  to  a free 
surface  (i.e.,  stress  or  pressure  equal  to  zero).  This  enables 
the  code  to  handle  a wide  variety  of  stress  wave  problems 

Each  layer  may  exhibit  different  material  properties  as 
characterized  by  a stress-specific  volume  relationship  accom- 
modating elastic,  plastic  and  compactible  ranges  as  shown  in 
Fig.  2.  This  relationship  is  represented  by  a piecewise  linear 
approximating  function  consisting  of  a number  (maximum  of  30) 
of  connected  straight  line  segments.  Provision  is  also  included 
to  allow  for  unloading  from  any  stress  state.  A maximum  of 
four  linear  unloading  paths  can  be  specified  for  various  ranges 
of  the  function. 

The  computer  code  can  handle  two  types  of  applied  loadings. 
First,  the  loading  upon  the  top  surface  of  the  material  may  be 
an  arbitrary  pressure- time  function  which  is  represented  in  the 
form  of  a stepwise  series  of  values  as  shown  in  Fig.  3,  Up  to 
75  steps  may  be  specified,  thereby  providing  a relatively  good 
fit  to  most  arbitrarily  desired  curves  (at  least  to  the  extent 
required  by  the  SWIMM  Code).  The  second  method  of  surface 
loading  is  supplied  by  the  specification  of  the  velocity  of  a 
second  moving  body  which  impacts  the  layered  media.  In  this 
manner  problems  such  as  a projectile  striking  a fixed  surface 
may  be  analyzed. 

The  wave  diagram  in  Fig.  4 presents  a qualitative  picture 
of  the  propagation  of  various  wave  fronts  in  the  layered  media. 
The  numbered  areas  indicate  distinct  zones  of  a uniform  state. 
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30  Segments  maximum 
4 Unloading  Curves 
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FIG.  2 STRESS-SPECIFIC  VOLUME  CURVE  FOR  A TYPICAL  LAYER 


FIG 


When  SWIMM  Code  is  run  on  the  IBM  7094  computer,  storage  capacity 
considerations  Limit  the  number  of  such  zones  which  can  be 
analyzed  to  approximately  1000. 


SHALLOW  Code  (Ref. 


A hydrodynamic  model  was  developed  for  predicting  the 
ground  shock  environment  from  surface  and  shallow-buried 
nuclear  explosions.  In  general  for  a shallow-buried  explosion 
at  early  times  the  direct  ground  shock  will  be  completely  con- 
tained as  indicated  in  Fig.  5.  A computer  code  was  written  to 
solve  the  one-dimensional  problem  of  determining  the  environ- 
ment at  the  shock  and  within  the  shock-engulfed  region.  The 
soil  can  have  a completely  arbitrary  Hugoniot  relation  (the 
loading  phase  of  the  pressure-volume  curve  of  Fig.  6.),  but  it 
is  restricted  to  expand  (unload)  along  the  same  curve.  This 
implies  that  the  soil  behaves  essentially  as  a nonlinear, 
elastic  media.  An  additional  assumption  is  that  the  material 
density  distribution  behind  the  shock  has  a prescribed  form.  The 
specific  form  used  for  the  model  is  given  in  Fig.  7.  The  incor- 
poration of  these  assumptions  into  the  code  allows  the  computation 
of  pressure,  density  and  particle  velocity  at  and  behind  the  shock 
front . 


This  approach  was  used  primarily  to  determine  media 
behavior  at  the  shock  front  and  has  been  applied  successfully 
for  this  purpose.  However,  while  the  predicted  gross  character- 
istics of  the  shock  engulfed  region  behind  the  front  may  be 
reasonable,  the  use  of  this  model  is  not  recommended  when  local 
response  behind  the  front  is  of  importance. 

At  a certain  time  after  burst,  the  direct  ground  shock 
will  strike  the  surface  and  initiate  the  propagation  of  a 
rarefaction  front  into  the  previously  shocked  region  as  illus- 
trated in  Fig.  8a.  The  rarefaction  front  will  eventually 
catch  up  with  the  one-dimensional  shock  wave  and  from  that 
point  on  will  increase  the  decay  rate  of  the  shock  (Fig.  8b). 
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Volume  Ratio  v/vq 

FIC.  6 HUCONIOT  AND  EXPANSION  RELATION  FOR 
A SOIL  PARTICLE 
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Shock  Front 


(a)  Direct  Ground  Shock  Before  Rarefaction  Catch-Up 


(b)  Direct  Ground  Shock  After  Rarefaction  Catch-Up 


FIG.  8 GROUND  SHOCK  STRUCTURE 
AT  LATER  TIMES 
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With  the  use  of  additional  simplifications,  a method  was  developed 
for  surface  bursts,  to  predict  the  characteristics  of  the  rare- 
faction front  and  of  the  direct  ground  shock  front  after  rare- 
faction catch  up.  Calculations  for  behavior  behind  the  fronts 
are  not  made.  The  code  also  has  the  capability  of  determining 
the  air-induced  shock  front  parameters  assuming  there  is  no 
interaction  with  the  direct  ground  shock. 


MULTIPLE  Code (Ref.  3; 


This  is  a computer  code  for  calculating  the  crater  size 
produced  by  multiple  underground  nuclear  explosions.  It  permits 
the  computation  of  total  earth  removal  for  two  explosions 
horizontally  displaced  as  a function  of  the  displacement 
distance,  the  depth  of  burst,  and  the  explosive  yield.  Figure  9 
portrays  the  three  zones  which  must  be  considered  - the  shock 
engulfed  region  of  (1)  the  first  burst,  (2)  the  second  burst 
and  (3)  their  interaction.  Also  included  in  the  computer  model 
is  the  effect  of  lag  time  between  the  initiation  of  the  first 
explosion  md  the  second.  Very  general  earth  properties  can 
be  handled  (the  same  as  for  SHALLOW  Code).  They  can  be  specified 
in  either  analytical  or  tabular  form. 


The  computer  program  calculates  the  emergent  one-dimensional 
shock  history  and  shock  engulfed  field  for  each  of  the  explosions 
in  a manner  similar  to  that  employed  in  the  SHALLOW  Code  for  the 
contained  shock  case.  In  the  interaction  zone,  the  shock  details 
are  determined  by  reducing  the  problem  to  one  dimension  by 
means  of  simplifying  assumptions.  When  one  of  the  shock  fronts 
first  reaches  the  surface,  the  total  upward  non-negative  mom- 
mentum  of  the  soil  is  computed  and  a resultant  crater  size 
(not  profile)  is  estimated.  The  results  indicate  that  the 
interaction  of  the  two  shock  waves  produces  an  enhancement  of 
the  upward  momentum  and  consequently  of  the  earth  removal 
capability. 
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4.  Stress  Wave  Propagation  in  Regions  of  Variable 
Modulus  of  Deformation  (Ref,  4)* 

A series  of  analytic  solutions  was  obtained  for  a class 
of  air-induced  ground  shock  problems.  These  exact  solutions 
deal  with  layered  media  in  which  the  modulus  of  deformation 
of  the  soil  increases  with  depth  in  a prescribed  manner.  The 
analysis  includes  both  a one-dimensional  nonsteady  solution 
and  a two-dimensional  quasi-steady  solution,  Fig.  10.  The 
latter  solution  is  limited  to  the  superseismic  regions  and 
yields  details  of  both  the  horizontal  and  vertical  motions 

The  equation  of  state  (See  Fig.  11)  which  is  used  to 
describe  the  material  behavior  is  a linear  and  reversible  one, 
i.e.,  a constant  modulus  of  deformation  at  any  one  depth.  The 
density  of  the  media  is  assumed  constant  throughout  the  entire 
region.  In  the  one-dimensional  nonsteady  treatment  no  additional 
materials  characteristics  need  be  specified;  however  for  the 
two-dimensional  quasi-steady  treatment  the  shear  capacity  of 
the  media  must  be  examined.  As  a preliminary  part  of  this 
latter  treatment  the  stress  along  the  initial  disturbance  front 
(the  dilatation  wave)  was  determined  for  both  an  elastic  and  a 
hydrodynamic  media.  There  was  little  difference  between  the 
two  cases  except  in  the  extreme  regions  where  the  onset  of  sub- 
seismic  wave  propagation  occurs.  Thus  the  hydrodynamic  model 
was  used  to  derive  the  complete  wave  propagation  solution  for 
the  two-dimensional  quasi-steady  problem. 

The  problem  is  further  limited  to  small  strains;  however 
this  condition  is  generally  satisfied  with  most  soils  and  rocks 
in  the  pressure  ranges  of  general  interest.  Thus  the  problem 
is  linear  and  the  principle  of  superposition  is  used  Exact 
basic  solutions  have  been  obtained  (Fig.  12  and  13)  for  a unit 
surface  or  blast  wave  load  (i.e.,  a simple  step  pulse  shown  in 


Reference  4 is  presented  as  the  Appendix  to  this  report . 
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(b)  Two-Dimensional  Quasi-Steady  Motion 

FIG  10  WAVE  PROPAGATION  GEOMETRY 


Specific  Volume 


flON-OF-STATE  RELATION 


L6 


FIG.  12  SOLUTION  FOR  FOURTH  POWER  MODULUS  VARIATION  AND  STEP  PULSE  LOAD,  ONE-DIMENSIONAL  CASE 
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FIG.  13  SOLUTION  FOR  FOURTH  POWER  MODULUS  VARIATION 
AND  STEP  PULSE  LOAD,  TWO-DIMENSIONAL  CASE 


Fig.  14a)  and  a requirement  on  the  modulus  of  deformation  var- 
iation with  depth  was  determined  for  both  the  one-dimensional 
and  two-dimensional  problems.  The  result  was  that  a relatively 
simple  fourth  power  variation  with  depth  (shown  quantitatively 
in  Fig.  15)  was  required.  This  can  be  fitted  to  most  actual 
site  data  quite  readily.  Simple  integral  equations  were  then 
developed  for  the  general  case  of  an  arbitrary  surface  load. 
Specific  solutions  are  also  given  for  the  more  common  case  of 
a simple  exponential  surface  load,  Fig.  14b. 

The  analyses  were  further  extended  to  include  the  more 
realistic  case  where  a variable  modulus  surface  layer  overlaps 
an  infinite  homogeneous  base  layer  requiring  the  use  of  a 
modified  molulus  form  as  indicated  in  Fig.  16.  This  configuration 
gives  rise  to  distinct  solution  zones  bounded  by  the  reverber- 
ation of  the  initial  disturbance  within  the  surface  or  variable 
modulus  layer.  The  solution  was  carried  out  to  include  approxi- 
mately four  solution  zones  and  programmed  for  a digital  computer. 
The  zones  for  the  one-dimensional  problem  are  depicted  in  Fig.  17. 
The  surface  load  is  defined  for  this  calculation  as  the  sum  of 
three  exponential  terms  (to  coincide  with  the  usual  air  blast 
loading  profile) . 

B . Solid  Models 
1.  DISCRETE  Code  (Ref.  5)* 

The  DISCRETE  computer  program  was  developed  to  obtain 
approximate  time-history  solutions  for  the  free-field  vertical 
component  of  soil  displacement,  velocity  and  stress  for  repre- 
sentative Titan  II  missile  sites  subjected  to  surface  air  blast 
loading.  It  is  coded  in  FORTRAN  IV  language  and  utilizes  the 
IBM  7094  computer  to  obtain  numerically  integrated  solutions 
for  a discrete  system  approximation  of  a one-dimensional  con- 
tinuous system. 

The  soil  model  employed  is  that  of  a one-dimensional 
layered  media  with  quite  general  piecewise  linear  stress-strain 

-'Reference  5 was  prepared  for  WES  and  issued  as  a separate  report 
on  this  contract. 
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(a)  Step  Pulse 


(b)  Exponential  Load 
FIG.  14  SURFACE  LOADINGS 
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Depth 


FIG.  15  PRESCRIBED  MODULUS  VARIATION  FORM 
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FIG.  16  MODIFIED  MODULUS  FORM 
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Zone  VI 


E DIAGRAM,  MODIFIED  MODULUS  FORM,  EXPONENTIAL  LOAD 


characteristics  specified  for  each  material  layer.  Each  of  the 
layers  can  possess  independent  loading  and  unloading/reloading 
stress-strain  characteristics  as  shown  in  Fig.  18.  There  is 
also  provision  for  specification  of  a maximum  admissible  tension 
for  each  layer. 

The  idealized  model  (Fig.  19)  coded  by  the  computer  program 
considers  a vertical  column  of  soil  that  is  subjected  to  a known 
loading  function  at  its  upper  end  and  connected  at  the  lower  end 
to  a linear  elastic  base  material  of  infinite  extent  with  known 
properties.  The  continuous  soil  column  is  idealized  as  a dis- 
crete element  system  in  which  each  element  represents  a soil 
layer  with  constant  properties  throughout.  The  assumption  of 
constant  properties  is  reasonably  accurate  if  a sufficient  number 
of  layers  are  selected  for  the  soil  model. 

Consideration  of  spatial  attenuation  of  stress  with  depth 
is  included  as  an  option  in  the  program.  Attenuation  effects 
are  introduced  by  altering  the  one-dimensional  soil  parameters 
automatically  within  the  program.  The  parameter  values  are 
altered  based  on  an  approximate  solution  to  the  Boussinesq 
problem  of  static  loading  on  the  free  surface  of  an  elastic 
half  space. 

The  response  of  the  discrete  element  system  is  obtained 
by  solution  of  the  differential  equations  of  motion  using  a 
Runge-Kutta  method  of  numerical  integration.  Output  from  the 
program  includes  the  transient  stress,  acceleration , velocity 
and  displacement  of  arbitrarily  selected  nodes.  An  option  is 
also  provided  for  determining  the  shock  spectrum  and/or  tape 
output  of  the  motion  and  stress  variables  for  given  time  incre- 
ments at  the  nodes.  This  latter  information  can  be  used  directly 
as  free-field  data  input  to  soil-structure  interaction  computer 
codes . 
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2.  SLAM  Code  (Ref.  6 and  7) 


The  SLAM  Code  (Stress  Waves  in  Layered  Arbitrary  Media) 
was  designed  primarily  to  investigate  the  geometrically  complex 
configurations  usually  encountered  at  missile  sites  of  interest. 
This  code  has  two  primary  flexibilities  built  into  it  which  are 
required  to  satisfactorily  attack  the  generally  complex  problems 
of  interest.  The  first  concerns  the  means  of  treating  the 
complex  geometric  boundaries  typically  encountered  at  sites  of 
interest,  and  the  second  concerns  the  ability  to  treat,  within 
the  same  problem,  differing  materials,  each  exhibiting  signifi- 
cantly differing  properties. 

Geometric  Complexities 

The  SLAM  Code  is  based  on  the  finite  element  approach  to 
the  stress  wave  problem  of  interest  herein,  as  opposed  to  the 
more  usual  finite  difference  techniques  that  have  been  applied 
to  this  class  of  problem.  The  finite  element  approach  can  be 
considered  as  merely  an  alternate  method  of  obtaining  an  equiv- 
alent finite  set  of  ordinary  differential  equations  which  approx- 
imate the  actual  governing  system  partial  differential  equations 
of  motion.  If,  in  particular,  a finite  element  mesh  is  chosen 
which  yields  a uniform  spacing  of  nodes  throughout  the  field, 
the  set  of  differential  equations  obtained  for  the  motion  of  the 
nodes  is  entirely  analogous  to  those  obtained  by  the  more  usual 
finite  difference  techniques. 

The  disadvantages  inherent  in  using  the  finite  element 
method  are  exactly  the  same  as  those  encountered  in  any  of  the 
finite  difference  or  discrete  element  methods  typically  used 
for  this  class  of  field  problems.  The  principal  drawback  of  any 
discrete  method  of  analysis  (that  may  be  noted  from  the  results 
obtained  by  the  various  methods)  lies  in  the  apparent  oscillations 
that  are  induced  into  the  computed  system  response.  These  oscil- 
lations come  about  simpl}'  from  the  fact  that  the  actual  system 
under  study  contains  an  infinite  set  of  system  frequencies  where- 


27 


as  the  approximating  set  contains  only  a finite  set  of  these 
frequencies . 

The  primary  advantage  of  the  finite  element  method,  however, 
over  the  finite  difference  formulation  lies  in  the  means  of 
treating  the  complex  geometries  typically  encountered.  For 
example,  a general  complex  site  geometry  is  shown  in  Fig  20. 

Since  the  site  boundaries  and  interfaces  are  generally  non- 
uniform,  the  use  of  finite  difference  formulations  is  made  quite 
difficult  by  the  relative  uniformity  of  mesh  point  spacing 
required  by  the  method.  The  use  of  the  finite  element  method 
however,  completely  eliminates  this  problem,  since  the  spacing 
of  nodes  throughout  the  free-field  is  completely  arbitrary  and 
can  be  made  to  satisfy  the  system  geometry.  This  is  indicated 
in  Fig.  21 . 

As  a further  example  of  the  flexibility  allowed  by  tnis 
method,  consider  the  problem  shown  in  Fig.  22  which  consists  of 
a three  layered  material  zone  which  passes  through  a fauit  that 
can  be  used  to  approximate  the  system  as  shown  in  Fig  23  It 
may  be  noted  how  simply  the  finite  element  mesh  can  be  matched 
to  the  complex  boundary  and  interface  conditions  of  this  problem. 

As  another  example  of  the  flexibility  of  the  method,  a 
finite  element  mesh  for  a uniform  half-space  problem  is  shown 
in  Fig.  24.  It  is  clear  of  course  that  for  this  simple  problem, 
a finite  difference  formulation  can  be  used.  However,  in  this 
problem  it  is  desired  to  have  a relatively  fine  mesh  at  the 
several  ground  ranges  of  interest  (shown  by  the  circles  in  Fig.  24) 
to  allow  for  sufficient  description  of  the  motion  history,  while 
to  obtain  a long  enough  duration  of  the  results,  the  mesh  bound- 
aries must  be  placed  sufficiently  far  so  that  they  will  not 
interfere  with  the  computed  results  during  the  times  of  interest. 
The  uniformity  of  mesh  spacing  required  by  the  finite  difference 
method  would  generally  require  an  excessively  large  number  of 
node  points  in  the  mesh  leading  to  a longer  computer  running  t im-  . 


Decay  of  Pressure 
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20  GENERAL  CONFIGURATION  OF  SITES  OF  INTEREST 


Ground  Range  of  Interest 


FIG.  22 


SITE  CONFIGURATION  FOR  PROBLEM 
EXHIBITING  A FAULT  ZONE 
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FIG.  24  TYPICAL  ELEMENT  MESH  FOR  UNIFORM  HALF  SPACE  PROBLEMS 


With  the  finite  element  method,  the  mesh  can  be  easily- 
expanded  away  from  the  zones  of  interest  allowing  for  a reduced 
number  of  node  points  requiring  consideration  and  leading  to  a 
shorter  computer  running  time.  These  zones  between  areas  of 
primary  interest  and  the  area  boundaries  act  primarily  as 
mechanisms  to  absorb  the  wave  without  introducing  reflected 
waves  from  any  more  closely  spaced  boundaries.  Several  methods 
for  treating  boundaries  to  eliminate  reflections  have  been  pro- 
posed, but  none  have  been  found  to  be  satisfactory,  particularly 
for  the  nonlinear  problem. 

As  can  be  noted  from  the  meshes  of  Fig.  23  and  24,  the 
element  meshes  considered  in  the  SLAM  Code  consist  of  triangular 
and  rectangular  elements.  The  fine  rectangular  meshes  are  used 
in  the  zones  of  primary  interest  where  computed  response  is 
desired.  Coarser  grids  are  used  in  zones  away  from  these  areas 
while  the  triangular  elements  are  used  to  merge  these  areas 

Missile  site  configurations  can  be  divided  into  two  classes 
of  problems,  (1)  those  made  up  of  uniform  horizontally  bedded 
one,  two  and  three  layered  problems,  and  (2)  unusual  configura- 
tions. For  the  first  class,  no  further  mesh  considerations  need 
be  mentioned.  In  fact,  the  same  mesh  shown  in  Fig.  24  for  the 
uniform  half-space  problem  can  be  used  to  analyze  two  and  three 
layered  problems  by  merely  changing  the  property  data  for  the 
various  zones  of  interest.  The  second  class  of  problems  indicates 
the  wide  versatility  of  the  method  over  other  techniques  A 
typical  mesh  for  one  of  these  problems  is  shown  in  Fig  23. 

It  should  also  be  mentioned  that  the  actual  mesh  formation 
for  specific  problems  is  highly  dependent  upon  the  material 
properties  of  the  system.  The  usual  criteria  (encountered  in 
finite  difference  methods  also)  that  must  be  followed  is  that 
the  time  for  propagation  of  the  stress  pulse  across  the  element 
must  be  small  compared  to  the  response  times  of  interest  in  the 
particular  problem,  at  least  in  the  zones  of  interest 


The  second  primary  area  of  flexibility  built  into  the  SLAM 
Code  concerns  the  freedom  of  choice  that  is  available  for  choosing 
material  constitutive  relations.  As  can  be  seen  from  some  of 
the  previous  discussion,  a typical  site  configuration  may  consist 
of  many  different  materials,  each  in  an  arbitrary  orientation 
to  each  other,  and  each  possessing  significantly  differing  material 
properties.  Thus  the  SLAM  Code  has  been  designed  to  allow  for 
an  arbitrary  number  of  zones  of  material,  each  having  the  capa- 
bility for  selecting  material  properties  from  an  arbitrary  list 
of  constitutive  relationships. 

The  catalogue  of  material  descriptions  for  which  numerical 
results  can  be  obtained  are  as  follows: 

• Isotropic  Elastic  Material 

This  material  property  is  the  simplest  to  treat,  of  course, 
and  does  have  application  to  various  material  problems  of  interest 
in  the  lower  overpressure  regimes. 

• Anisotropic  Elastic  Materials 

This  capability  was  added  to  the  list  to  allow  for  inves- 
tigation of  various  rock  foundation  problems  which,  due  to  the 
method  of  formation  of  the  material,  have  differing  elastic 


properties  in  the  horizontal  and  vertical  directions.  These 
bedding  problems  can  be  easily  treated  in  this  manner  in  for- 
tunately, for  most  rock  material  of  interest,  little  property 
data  is  available  with  which  to  determine  the  necessary  constants. 
(Unfortunately,  this  problem  is  encountered  for  many  real  soil/ 
rock  materials . ) 


• Viscoelastic  Materials 

To  treat  materials  which  exhibit  significant  damping  char- 
acteristics (i.e.  clay  soils),  a simple  viscoelastic  constitutive 
model  has  been  incorporated  into  the  system. 
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• Compressible  Fluid 

Since  the  SLAM  Code  has  been  used  to  investigate  configura- 
tions in  which  stress  wave  propagation  from  solid  to  water  and 
vice  versa  was  of  interest,  a linear  compressible  fluid  model 
was  incorporated. 

• Elastic-Mises  Plastic,  Arbitrary  Strain 
Hardening  Material 

To  account  for  nonlinear  material  behavior,  one  of  the  con- 
stitutive relations  incorporated  into  the  code  is  an  elastic- 
plastic  material  in  which  plastic  flow  is  governed  by  the  Prandt- 
Reuss  relations  and  the  initial  yield  is  based  upon  the  Mises 
yield  criteria.  As  seen  in  Fig.  26,  the  yield  surface  considered 
is  a cylindrical  surface  in  principal  stress  space,  with  an  axis 
(Aline  of  Fig.  26a)  that  makes  equal  angles  with  the  principal 
stress  axes.  When  the  stress  state  reaches  a point  on  the  yield 
surface,  initial  plastic  flow  can  occur.  This  point  corresponds 
to  the  point,  a , on  the  effective  stress-strain  relation  of 
Fig.  26b.  As  further  straining  occurs  the  material  becomes 
stronger  (strain  hardens)  and  the  diameter  of  the  yield  surface 
increases.  For  unloading,  recovery  of  the  material  takes  place 
elastically.  The  effective  stress-strain  curve  of  Fig.  26b  can 
be  obtained  from  the  results  of  an  unconfined  compression  test. 

In  the  SLAM  Code,  any  arbitrary  strain  hardening  law  can  be  used 
as  data  input  to  the  code  by  considering  it  in  the  form  of  a 
series  of  points  along  the  effective  stress-strain  curve  of  the 
material . 

This  plastic  behavior  model  was  initially  developed  for 
metal  materials  in  which  two  primary  characteristics  are 
apparent:  (1)  the  material  behaves  similarly  in  both  tension 

and  compression,  and  (2)  no  plastic  volume  change  occurs  'hydro- 
static stress  states,  states  on  thevYline,  produce  only  elastic 
volume  changes).  Some  rock  materials  at  the  lower  overpressure 
ranges  do  exhibit  this  behavior  under  compressive  loadings 
However,  under  tensile  loadings,  they  exhibit  a much  different 
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(a)  Yield  Surface  for  Mises  Yield  Condition 


(b)  Strain  Hardening  Characteristics 


FIG  26  ELASTIC- PLASTIC  BEHAVIOR  OF  MISES  MATERIAL 


behavior  so  that  plastic  flow  occurs  even  under  hydrostatic  tensile 
loadings.  This  behavior  can  be  incorporated  into  the  model  by 
using  a spherical  cap  on  the  tension  side  of  the  yield  surface 
(Fig.  26a). 

• Simple  Coulomb-Mohr  Plastic  Flow  Theory 

For  some  rock  and  most  soil  materials  the  Mises  theory  is 
inapplicable  since  the  materials  exhibit  a significant  increase 
in  strength  with  hydrostatic  compressive  pressures.  To  include 
this  effect  into  the  material  model,  a plastic  yield  law  was 
developed  which  incorporated  this  effect.  This  yield  criterion, 
plotted  in  the  same  principal  stress  space,  is  a conical  surface 
with  an  axis  again  ir.  theAline,  as  seen  in  Fig.  27.  For  cohesion- 
less materials  (sands) , the  apex  of  the  cone  is  located  at  the 
(0,0,0)  point  in  stress  space  while  the  opening  angle  of  the 
cone  is  related  to  the  coefficient  of  friction  of  the  soil /rock 
which  is  typically  obtained  from  the  triaxial  compression  test . 

For  cohesive  materials,  the  location  of  the  apex  along  the Aline 
is  related  to  both  the  coefficient  of  friction  and  the  cohesion 
obtained  from  the  triaxial  test. 

• Modified  Coulomb-Mohr  Plastic  Flow  Model 

On  a current  BSD  study  being  conducted  at  IITRI,  improved 
models  of  plastic  yield  surfaces  are  being  investigated  which 
attempt  to  more  closely  link  the  behavior  of  real  soils  or  rocks 
to  the  analytic  constitutive  formulations.  An  example  of  this 
is  shown  in  Fig.  28  in  which  the  same  Coulomb-Mohr  yield  surface 
is  used  as  before.  However,  as  may  be  seen  in  the  figure,  the 
conical  yield  surface  is  capped  by  spherical  surfaces  at  be  th 
the  compression  side  and  the  tension  side.  This  enables  two 
relatively  well  known  properties  of  soils  to  be  simulated 

First,  the  result  of  a series  of  triaxial  tests  on  a typ- 
ical cohesive  soil  is  shown  in  Fig.  29.  The  linear  Mohr  envelope 
applies  at  relatively  low  values  of  confining  stress  ievels. 

But,  from  the  results  of  an  unconfined  triaxial  test  (no  lateral 
restraining  pressures),  the  cohesive  strength  C^,  is  usual lv 
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much  lower  than  the  apparent  cohesion,  C,  obtained  by  extrapolating 
the  envelope  to  the  zero  stress  line.  By  placing  the  cap  on  the 
conical  yield  surface  near  the  apex  of  the  cone,  this  effect  can 
be  relatively  easily  simulated. 

On  the  other  hand,  by  using  the  cap  at  the  compressive 
side  of  the  cone  one  can  then  simulate  the  nonlinear  behavior 
of  soils  noted  in  the  hydrostatic  and  confined  (consolidation) 
tests.  A typical  example  of  this  plastic  behavior  during  confined 
compression  is  shown  in  Fig.  30.  The  use  of  the  simple  Coulomb- 
Mohr  theory  alone  would  not  allow  the  simulation  of  this  complex 
behavior.  By  allowing  motion  of  the  compressive  cap  with  plastic 
loading  and  unloading,  the  complex  stress  history  of  Fig  30  can 
be  approximated. 

3.  Deviatoric  Effects  in  High- Intensity  Stress  Waves  jRef 8)_ 

A study  was  made  of  the  one-dimensional  propagation  of 
strong  spherical  shock  waves  in  solids.  The  phenomena  associated 
with  cratering  and  direct  ground  shock  due  to  nuclear  detonations 
were  of  prime  concern. 

The  problem  posed  and  analyzed  is  that  of  a media  containing 
a spherical  cavity  which  is  initially  loaded  by  extremely  high 
pressure  on  its  surface.  A shock  wave  is  created  in  the  surround- 
ing solid  forming  successive  zones  of  primarily  fluid,  plastic 
and  elastic  behavior  as  the  wave  propagates  outward  and  decays 
(Fig.  31).  The  fluid  and  elastic  zones  have  been  studied  in 
great  detail  whereas  the  plastic  zone  has  been  given  relatively 
little  attention.  The  present  analysis  bridges  the  gap  between 
the  fluid  and  elastic  solutions  by  consideration  of  plastic 
behavior.  In  fact  it  affects  a continuous  transition  from  one 
zone  to  another  by  using  a material  description  which  accommo- 
dates all  three  types  of  behavior  together.  In  the  formulation, 
the  plastic  (deviatoric)  effects  remain  the  same  order  of  magni- 
tude regardless  of  stress  level;  hence  at  extremely  high  stress 
levels  a fluid-like  behavior  dominates.  The  inclusion  of 
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FIG.  30  STRESS- STRAIN  RELATION  FOR  OTTAWA  SAND 


ZONES  OF  BASICALLY  DIFFERENT  BEHAVIOR 


deviatoric  effects  required  the  development  of  a rational  descrip- 
tion of  real  media  applicable  to  any  three-dimensional  state  of 
stress.  This  was  done  by  extending  extablished  plasticity  theo- 
ries to  include  compressibility  effects,  and  by  modification  of 
present  hardening  concepts  to  conform  more  nearly  tc  the  at tual 
physics  of  material  behavior. 

This  analysis  has  demonstrated  that  the  deviatoric  effects 
can  aLter  the  response  in  the  very  close-in  region.  In  essence 
the  presence  of  deviatoric  stress  permits  the  radial  transmission 
cf  a compressive  stress,  despite  the  fact  that  the  material  ex- 
periences a hydrostatic  tension. 

Techniques  for  solving  the  governing  equation^  hive  been 
outlined.  A combination  of  numerical  and  analytical  procedures 
was  indicated  to  integrate  the  differential  equations.  At  'he 
shock  front  the  Rankine-Hugoniot  relations  are  used  tcg-Ltiei 
with  a progressive  wave  solution  to  analytically  extend  the  so- 
lution into  a relatively  smooth  region.  The  remainder  ci  tne 
solution  (behind  the  shock  front)  can  then  be  determined  rum  t 
icaliy . 
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III.  SOIL- STRUCTURE  INTERACTION 


A.  VIM  Code  I (Ref.  9) 

This  program  was  written  to  determine  the  vertical  rigid- 
body  motion  imparted  to  the  contents  of  a silo  subjected  to  air 
blast  surface  overpressure.  In  the  analysis  the  silo  is  idealized 
as  a stepped,  rigid,  circular  cylinder  containing  two  constant- 
diameter  sections  to  simulate  the  Titan  II  Silo  configuration 
as  shown  in  Fig.  32. 

The  forces  acting  on  the  silo  which  were  considered  in  the 
^rtical  response  analysis  are  the  vertical  cap  force  F't),  skin 
friction  , vertical  loads  at  the  base  Fg,  and  the  soil -structure 
interaction  forces  Fp  at  the  step  in  the  cylinder 

The  time  varying  cap  load  was  taken  to  be  simply  the  speci- 
fied air  blast  overpressure  a times  the  cap  area.  Skin  fri.tion 
forces  were  computed  using  the  Coulomb  friction  model  which  utilizes 
information  on  free-field  soil  stresses  and  velocities  at  selected 
node  points  along  the  length  of  the  silo.  The  silo  base-soil 
interaction  vertical  resultant  force  Fg  was  calculated  on  the 
basis  of  a commonly  used  soil-structure  interaction  model  see 
Fig.  32).  This  model  stipulates  that  the  total  force  is  composed 
of  the  sum  of  three  forces  due  to  (1)  the  free-field  stress, 

(2)  the  relative  displacements  of  the  free-field  soil  and  the 
structure,  and  (3)  the  relative  values  of  the  free-field  soil 
velocity  and  the  velocity  of  the  structure.  The  constant  pro- 
portionality factor  kg,  (spring  constant)  for  the  relative  dis- 
placements is  assumed  to  be  equal  to  a fraction  of  the  spring 
constant  corresponding  to  the  elasticity  solution  for  the  dis- 
placement of  a rigid  circular  die  pressed  into  the  plane  boundary 
of  an  elastic  half  space.  The  calculation  of  the  base  damping 
constant  Sg,  associated  with  the  interaction  forces  produ^ ed  by 
the  relative  velocity  between  the  silo  and  the  free-field  soil 
is  also  based  on  an  elasticity  solution.  It  is  taken  to  be  a 
fraction  of  the  normal  stress  produced  at  the  end  of  a bat  which 
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is  struck  by  a rigid  body  with  a prescribed  velocity  relative 
to  the  bar. 

The  interaction  model  used  to  compute  the  vertical  force 
Fp  acting  on  the  underside  of  the  step  in  the  silo  is  similar 
to  that  for  the  silo  base,  with  the  exception  that  the  spring 
and  damping  constants  are  obtained  differently.  The  spring  con- 
stant kp  is  based  on  the  average  displacement  of  a uniform  ver- 
tically loaded  annular  region  acting  on  the  plane  boundary  of 
an  elastic  half  space  shown  in  Fig.  33.  The  damping  constant 
Sp  is  obtained  by  assuming  that  the  ratio  of  damping  to  spring 
constant  is  the  same  at  the  step  point  in  the  silo  and  at  the 
silo  base. 

The  equation  of  veritcal  motion  for  the  rigid  silo  is 
integrated  by  the  Runge-Kutta  method  of  integration  attributed 
to  Gill.  Output  from  the  program  includes  the  verti. ai  displace- 
ment, velocity  and  acceleration  of  the  silo;  the  interaction 
stresses  and  forces  at  the  base  and  step  positions;  and  the 
skin  friction  stresses  along  the  length  of  the  silo 

B.  VIM  Code  II  (Ref.  10) 

This  code  treats  the  vertical  interaction  of  a silo  with 
a stress  wave  passing  through  the  soil.  It  is  a generalization 
of  VIM  Code  I and  was  written  specifically  to  handle  the  single 
silo  structure  or  the  "articulated"  silo  structure,  as  tvpified 
by  the  LER/Launch  Tube  construction  employed  with  the  Minuteman 
Launch  Facility. 

The  vertical  flexibility  of  the  silo  structure  is  accounted 
for  by  considering  it  as  a set  of  discrete  mass  points  inter- 
connected by  stiff  springs,  the  stiffnesses  of  which  are  determined 
by  the  cross-sectional  and  elastic  properties  of  the  silo,  The 
idealization  is  illustrated  in  Fig.  34. 

The  interaction  model  includes  effects  of  silo  skin  friction 
S^..  The  skin  friction  model  used  is  based  on  the  Coulomb  soil 
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shear  failure  theory.  The  frictional  forces  on  the  silo  are 
computed  from  free-field  inputs  which  are  described  at  the  mid- 
point of  soil  "layers"  in  terms  of  both  vertical  velocity  and 
horizontal  stress  pulses,  the  layers  being  chosen  so  that  the 
free-field  pulses  can  be  considered  as  fairly  uniform  over  ^ach 
layer.  The  free-field  data  can  either  be  read  in  or  determined 
from  idealized  motion  pulse  shapes  which  are  obtained  by  fitting 
pulses  to  shock  spectra. 

The  air  overpressure  load  FA(t)  on  the  cap  of  the  silo  is 
computed  as  the  value  of  the  overpressure  times  the  cap  area, 
neglecting  early  time  cap  engulfment  effects.  The  loading  (t) 
on  the  base  of  the  silo  is  determined  by  use  of  a soil-structure 
interaction  law  which  includes  the  soil  foundation  modulus  and 
damping  coefficient  which  are  prescribed. 

In  the  case  of  articulated  silos,  the  base  loads  Fg^'.'t) 
acting  on  the  upper  section  of  the  silo  are  computed  by  assuming 
that  a wedge  of  soil  fails  under  its  foundation.  These  leads 
are  incorporated  into  the  code  in  the  form  of  a spr ing-dashput 
model,  the  parameters  of  which  must  be  specified. 

The  code  produces  output  which  includes  the  displacement, 
velocity  and  acceleration  of  each  silo  mass  point,  as  well  as 
the  associated  soil  displacement  and  velocity  at  each  specified 
value  of  time.  The  computer  program  is  written  in  FORTRAN  IV 
language  for  the  IBM  7094  computer. 

C.  RING  Code  (Ref.  11) 

This  code  was  written  to  predict  the  response  of  buried 
flexible  tunnel  liners  encased  in  crushable  foam  materials  and 
subjected  to  the  engulfment  of  a plane  compression  wave  moving 
normal  to  the  axis  of  the  tunnel.  This  is  demonstrated  in  Fig.  35. 

The  soil-structure  interaction  model  on  which  the  code  is 
’i'.  d uncouples  the  free-field  response  from  the  structural 
r ->ponse . The  code  accepts  as  input  free-field  stress  and/or 
motion  history  data  directly  in  the  form  of  tabulated  values  at 
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various  times,  or  indirectly  as  parameters  of  an  ideal  shock 
"fitted"  to  actual  spectra  data.  The  soil  is  characterized  by 
a power-law  stress  strain  curve  suggested  by  experiments  and 
provides  for  a hysteresis  effect  during  unloading,  see  Fig.  36. 

In  the  interaction  model  the  free-field  stresses  and  those  gen- 
erated by  the  relative  displacement  and  velocity  of  the  soil  and 
structure  constitute  the  applied  load  to  the  outer  (foam)  structure. 
Only  the  normal  stress  components  are  considered  to  be  acting  on 
the  structure  - frictional  forces  are  ignored.  Stress  wave  effects 
through  the  foam  are  also  neglected.  The  most  general  type  of 
foam  that  can  be  presently  accommodated  by  the  code  is  that  which 
has  a triiinear  stress-strain  relation  (Fig.  37)  corresponding 
to  linear  elastic,  crushing  and  hardening  phases  of  behavior, 
with  provision  for  nonrecoverable  plastic  straining  upon  unloading. 
The  tunnel  liner  is  treated  as  a linear  elastic  ring  which  is 
analyzed  by  modal  techniques.  The  prediction  of  the  dynamic 
response  of  the  liner  involves  the  numerical  integration  of  the 
system  of  differential  equations  of  motion  which  is  performed 
by  use  of  the  Runge-Kutta-Gill  method. 

The  computer  program  is  written  in  the  FORTRAN  IV  language 
for  the  IBM  7094  digital  computer,  using  the  standard  1BSYS  monitor 
(version  13)  system.  The  output  information  that  can  be  obtained 
from  the  program  includes  the  time  history  as  well  as  the  maxima 
of  several  response  variables  at  selected  points  in  the  soil, 
foam  and  tunnel. 

D.  SIM  Code  (Ref.  12  and  13) 

The  SIM  Code  was  written  to  determine  the  motions  of  a 
missile  isolated  from  a hardened  silo  structure  which  is  sub- 
jected to  ground  shock  developed  by  a nuclear  detonation  as 
depicted  in  Fig.  38.  The  analysis  developed  for  the  system 
considers  beam  motions  of  both  the  missile  and  silo.  Although 
the  analysis  is  primarily  concerned  with  horizontal  response 
produced  by  horizontal  loads,  the  effect  of  rigid  body  vertical 
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FIG.  37  ASSUMED  STRESS-STRAIN  CURVE  FOR  FOAM 
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motion  of  the  silo  on  this  response  can  be  determined  from  the 
program  also. 

The  soil-structure  interaction  model  employed  separates 
the  soil  and  structural  response  problems.  The  model  is  the 
same  as  in  the  RING  Code.  The  loads  imposed  on  the  silo  are 
those  assc  iated  with  the  free-field  stress  and  the  relative 
displacements  and  velocities  of  the  soil  and  structure  The 
free-field  motion  histories  with  depth  can  either  be  r -ad  in 
or  computed  internally  from  specified  shock  spectra  parameters. 

The  free-field  stress  is  calculated  from  the  motion  data  by 
assuming  that  the  soil  at  each  depth  is  in  a state  of  plane 
strain . 

The  silo  and  internal  components  are  idealized  as  a number 
of  interconnected  elastic  flexible  members  (up  to  a maximum  of 
five)  of  which  the  first  is  the  primary  or  silo  structure  with 
the  other  housed  within.  The  idealization  is  indicated  in  Fig  39. 
Each  member  of  this  system  is  assumed  to  respond  as  a flexible 
beam  structure  including  the  effects  (if  desired)  of  shear  defor- 
mation and  rotary  inertia  along  each  member.  The  beams  are 
analyzed  as  lumped  mass  systems  (Fig.  40)  with  up  to  60  degrees 
of  freedom.  It  is  assumed  that  in  the  plane  of  the  silo  cross 
section,  the  silo  is  rigid,  that  is,  ring  bending  eftevts  have 
b-en  ignored  as  illustrated  in  Fig.  41. 

The  computer  program  written  for  the  system  of  beams  allows 
for  general  types  of  connections  between  any  beams  at  any  voca- 
tions along  the  beams  The  types  of  interbeam  connections  cur- 
rently allowed  for  in  the  computation  scheme  have  been  geared 
for  a particular  application,  although  any  type  of  interconnection 
may  be  considered.  For  any  one  connection  the  general  model 
us  d is  representative  of  a wide  class  of  practical  connectors 
a J onsists  of  a nonlinear  spring  element,  a linear  dashpot 
aid  a number  (up  to  five)  of  Maxwell  elements  in  paraiLel  «Fig.  42a). 
i addition  to  these  force  transmission  connectors,  arbitrary 
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FIG.  39  IDEALIZATION  OF  SILO  SYSTEM 
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(a)  Force  Transmission  Connector 
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(b)  Moment  Transmission  Connector 


FIG.  42  TYPES  OF  INTERBEAM  CONNECTIONS  ALLOWED  IN  CURRENT  PROGRAM 
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moment  connections  are  provided  for  in  the  form  of  simple  torsion 
spring-dashpots  (Fig.  42b). 

The  computer  program  was  written  in  FORTRAN  TV  language 
for  the  IBM  7094  digital  computer.  The  input  required  to  the 
code  can  be  generally  grouped  into  three  categories:  beam  data, 

connection  data  and  ground  motion  data.  The  typical  output  from 
the  code  consists  of  the  displacement,  velocity  and  acceleration 
of  each  beam  at  various  points  along  the  beam,  applied  loads  and 
moments  to  the  members,  and  the  shear  and  moment  developed  in 
each  beam  along  its  length.  A further  output  option  is  provided 
to  the  user  in  that  this  data  can  be  printed  out  at  various  in- 
tervals through  the  course  of  the  system's  motion,  or  only  maximum 
values  can  be  printed  at  the  end  of  the  run. 
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x,  INTRODUCTION 


The  study  of  the  propagation  of  stress  waves  in  soils  and 
rocks  is  vital  to  the  design  of  buried  protective  structures 
against  nuclear  attack.  The  in  situ  soil  properties  are  generally 
poorly  defined,  however,  some  general  characteristics  have  been 
established.  Wilson  and  Sibley  (Ref.  14)  have  developed  a pro- 
cedure for  selecting  the  constrained  modulus  and  the  strain 
recovery  for  typical  surface  soil  conditions.  The  constrained 
modulus  can  increase  rapidly  with  depth  in  several  hundred  feet 
and  change  by  several  orders  of  magnitude.  The  strain  recovery, 
on  the  other  hand,  varies  from  approximately  60  percent  at  the 
surface  to  100  percent  at  depths  of  the  order  of  100  ft.  The 
density  of  the  soil  is  relatively  constant,  increasing  slightly 
with  increasing  depth.  Figure  43  presents  some  data  for  two 
actual  construction  sites. 

The  following  analysis  was  performed  in  order  to  obtain 
an  estimate  of  the  effect  that  the  region  of  increasing  modulus 
of  deformation  has  upon  the  variables  of  the  transient  stress 
field  which  results  from  the  application  of  a suddenly  applied 
load,  such  as  is  experienced  by  the  soil  during  a nuclear  attack. 
The  following  model  was  selected  because  its  specific  properties 
permits  one  to  obtain  simple  analytical  expressions  and  thus 
rapidly  provide  approximate,  but  reasonable  estimates  of  the 
soil  motions. 

This  paper  treats  both  one-dimensional  r.c-r.s:  s-  -dy  wave 
propagation  into  a region  of  changing  modulus  of  deformation 
and  a corresponding  two-dimensional  quasi-steady  problem.  This 
paper  will  summarize  the  procedures  used  and  the  results  obtained 
recently  by  the  author.  A more  detailed  discussion  is  given  in 
two  recent  papers,  Ref.'  15  and  16. 


II. 


ONE -DIMENSIONAL  NONSTEADY  STRESS  WAVE  PROPAGATION 


This  section  deals  with  one-dimensional  propagation  of 
stress  waves  into  a region  of  changing  modulus  of  deformation. 

The  constrained  modulus  of  deformation  of  the  soil,  M,  is  assumed 
to  vary  with  depth,  x,  (measured  from  the  surface),  in  any  manner. 
More  specific  variations  of  the  modulus  will  be  treated  in  sub- 
sequent sections  of  this  paper.  Furthermore,  the  density  of  the 
soil  is  assumed  to  be  constant  at  a value  pQ.  This  latter  assump- 
tion is  not  essential  to  the  development  of  some  of  the  results 
which  follow,  however,  it  does  permit  considerable  simplification. 
A linear  and  reversible  equation  of  state  or  stress-strain  law 
was  selected. 

In  general,  we  will  consider  a semi-infinite  half  space 
which  is  subjected  to  impulsive  loads  (stress)  at  the  surface 
(x=o) . Disturbances  will  propagate  down  into  the  medium  and 
interact  with  the  moduli  gradients.  The  assumption  of  a linear 
and  reversible  stress-strain  law,  together  with  a seismic  velocity 
which  is  dependent  upon  position,  leads  to  characteristics  or 
disturbance  paths  in  the  wave  diagram  which  are  self  similar. 
Furthermore,  the  initial  disturbance  or  shock  front  is  uncoupled 
from  the  interior  of  the  wave  as  well  as  from  subsequent  distur- 
bances at  the  boundary. 

The  interaction  of  a stress  wave  with  an  interface  across 
which  the  modulus  changes  from  M to  M+AM  will  result  in  the 
transmission  and  reflection  of  two  stress  waves.  The  incident 
and  reflected  wave  will  propagate  at  a velocity  c,  whereas  the 
transmitted  wave  will  propagate  at  a velocity  c+ac  corresponding 
to  the  modulus  M+AM. 

The  characteristic  equations  for  each  of  the  waves  can 
be  solved  and  one  obtains  the  following  results. 
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AM 
M ’ 


(AM  <<  M) 


(2) 


where  L r and  Au  are  the  changes  in  stress  and  particle 
velocity  across  the  interface. 

In  the  limit  as  aM  — 0,  Eq . (1)  becomes 


du  ^ o C u 

which  when  integrated,  yields 
cu  = Constant  = uQ 

where  a and  u are  values  from  some  reference  state  S , say 
0 0 “ 
the  initial  value  at  the  surface  (x=o) . The  hodograph  diagram, 

Fig.  44  illustrates  the  locus  of  these  shock  front . states . 

Integrating  directly  the  limiting  forms  of  Eq.  (2)  yields 


where  c and  M refer  to  state  S . 
o o o 

It  will  be  of  interest  to  compare  the  results  of  Eq , (3) 
with  the  case  of  a modulus  discontinuity  of  finite  magnitude,  viz 
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which  in  the  limit  for  a rigid  boundary  (M^  — ► ) yields  a 

reflection  coefficient  (for  stress)  of  or  / a =2.  The 

case  of  a continuously  varying  modulus  from  Mq  to  M1  yields 
much  larger  reflection  coefficients  for  both  stress  and  particle 
velocities.  If  / Mq  = 100,  which  is  not  uncommon  in  geologic 
materials  then  /oq  = 3.16  and  in  the  limit  as  , 

oo  . This  comparison  is  illustrated  in  Fig.  44. 

The  results  of  the  previous  paragraphs  permit  us  to  compute 
the  initial  state  behind  the  shock  front  as  it  propagates  in  the 
nonhomogeneous  medium.  This  information  is  of  value  since  the 
maximum  stresses  and  velocities  generally  occur  at  the  shock 
front  for  many  practical  problems.  However,  for  this  same  class 
of  problem  the  details  behind  the  front  are  also  of  considerable 
interest . 

In  numerically  examining  the  stress  field  behind  a shock 
wave  moving  into  a region  of  increasing  modulus  and  generated 
by  a step  surface  load,  it  appeared  that  the  particle  velocity 
was  relatively  independent  of  time.  This  suggested  that  perhaps 
there  exists  a variation  of  modulus  M(x)  for  which  the  particle 
velocity  is  only  a function  of  time.  We  will  seek  to  determine 
such  a modulus  variation.  Since  the  problem  is  linear  and  the 
principle  of  superposition  is  valid  it  is  only  necessary  to 
determine  the  response  of  the  medium  to  a step  load  (stress) 
applied  to  the  surface.  The  response  of  the  medium  to  a more 
general  surface  load  can  then  be  obtained  by  suitable  integra- 
tion processes. 


The 

field  equations  are: 

+ i-  o 

(4) 

r t £ 

Po 
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(5) 
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FIG.  44  INTERACTION  OF  STRESS  WAVE  WITH  MODULUS  GRADIENT 
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If  we  assume  that: 


r = a (x) 

and  substitute  this  form  into  Eq . (4)  and  (5)  we  obtain 


r = a + ax 

o — 

and 

u - u + — 

° Po 

where  a is  a constant. 

The  modulus  variation  must  be 

^ - <1  + f>  * U ± O 

where  b is  a constant.  Time  t is  made  dimensionless  by  defining 
t = t cQ/b.  This  basis  solution  is  illustrated  in  Fig.  45  and  46. 
The  inverse  solution  will  be  used  for  treating  reflected  wave 
systems  where  the  stress  waves  propagate  into  a region  of  decreas- 
ing modulus.  The  displacements  can  be  computed  by  directly 
integrating  the  velocity  with  respect  to  time.  The  shock  arrival 
time,  i , is  given  by 

- t+t 

We  will  now  consider  the  more  general  boundary  condition 
where  the  stress  at  the  surface  is  equal  to  some  function  G (■)• 
The  stress  varies  as 

ff-LL.  ll.  ~ (i  + n c(  t-t'> 
ao 

that  is  the  stress  is  magnified  by  the  factor  (1+0  ana  has  the 
same  normalized  shape  as  the  applied  load  at  the  boundary  but 
delayed  by  the  shock  arrival  time  t' . 

The  general  form  for  the  particle  velocity  is  not  quite 
as  simple.  We  can  write  the  velocity  in  the  following  form. 
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45  SOLUTION  FOR  FOURTH  POWER  MODULUS  VARIATION 


r 


— * ( 1 - T ' ) - ( T-T ' ) 

uo 


The  first  term  is  the  contribution  due  to  the  intensity  of  the 
surface  load  (or  initial  value)  and  the  second  term  is  the  time 
decay  contribution  which  is  delayed  by  the  shock  arrival  time  t' 
In  general,  then  the  velocity  will  consist  of  two  parts,  viz. 


u 

= (1  * T ' ) C(  T-T')  ^ / C(  -V- 


: ) d't 


field 

loads 

early 


The  dynamic  surface  loads  which  are  of  interest  to  the 
of  protective  construction  are  suddenly  applied  decaying 
These  loads  can  be  approximated  very  well,  at  least  for 
times,  by  a simple  exponential  form,  viz. 


^ / \ Ul 

G ( t ) e 

where  n is  a decay  factor.  The  previously  indicated  integra- 
tions can  be  carried  out  to  yield  the  following  results 
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Some  typical  displacement  results  are  shown  in  Fig.  47.  The 

response  of  the  soil  exhibits  a definite  bounce  characteristic. 

The  effect  of  the  decaying  surface  wave  decreases  with  depth 

since  the  response  of  the  soil  at  depth  increases.  For  a step 

pulse  the  soil  reaches  a maximum  displacement  s which  falls  off 

m 

as  the  inverse  square  of  the  depth  variable  , 


The  fourth  power  fit  of  the  modulus  variation  appears  to 
be  a reasonable  approximation  for  surface  and  near  surface  soil 
conditions.  However,  at  greater  depths  the  modulus  generally 
approaches  a constant  value.  The  results  of  the  previous  para- 
graphs can  be  used  to  construct  relatively  simple  analytical 
relations  to  describe  the  response  of  a media  of  this  more  com- 
plicated modulus  distribution.  To  illustrate  this  we  have 
selected  a modified  modulus  form  which  consist  of  a fourth- 
power  increase  to  a depth  t ^ and  then  a region  of  constant 
modulus,  . This  is  illustrated  in  Fig.  48.  We  will  derive 
the  response  of  the  medium  to  a step  load  at  the  surface  of 
intensity  r . Again  the  response  of  the  medium  to  a more  gen- 
eral input  G (•;-)  can  be  obtained  by  suitable  integrations. 

The  propagation  of  the  shock  front  will  occur  as  in  the 
basic  solution  and  result  in  the  generation  of  a reflected  and 
transmitted  disturbance  at  f = C j . This  disturbance  will 
reverberate  between  f =0  and  \ = [ i and  define  a number  of 
distinct  solution  zones.  These  are  illustrated  and  defined 
in  Fig.  49.  It  should  be  nested  that  the  zones  in  the  constant 
modulus  portion  of  the  medium  are  regions  of  simple  waves.  This 
fact  results  in  a simple  relationship  between  the  stress  and 
thj  particle  velocity. 
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FIG.  49  WAVE  DIAGRAM,  MODIFIED  MODULUS  FORM 


The  solution  has  been  constructed  for  the  first  four  zones. 
Typical  displacements  are  presented  in  Fig.  50. 

It  can  be  shown  using  the  requirement  of  conservation  of 
momentum  that  the  long  time  behavior  of  a semi-infinite  region 
of  constant  modulus  , covered  by  a relatively  thin  region  of 
varying  properties  is  governed  by  the  properties  of  the  constant 
region.  It  is  of  interest  then  to  examine  the  response  of  the 
medium  in  the  absence  of  the  surface  layer.  As  before  we  will 
consider  the  simple  case  of  a suddenly  applied  constant  load  of 
intensity  cQ  . A stress  wave  of  intensity  a will  propagate 

into  the  medium  at  a velocity  c,  and  accelerate  the  medium  to 

* , L 
a velocity  u , where 


P ci 
*o  1 

The  displacement  s at  the  surface  of  this  uniform  region  will 


3 _ -) 

T~  ' 2 1 + 

O 1 


At  late  times,  when  the  disturbances  reverberating  in  the 

surface  region  decay  the  entire  region  will  approach  a static 

• ^ 

condition  (compressed  by  a stress  c ) moving  at  a velocity  u 
The  displacement  of  particles  in  the  variable  modulus  region 
(0  < C,  ^ Cj  ) will  approach  that  given  by  Eq.  (6)  plus  an  incre- 
ment of  displacement  As'  which  is  due  to  the  compression  of  the 
entire  media  to  a stress  cr  . This  increment  of  displacement 
can  be  obtained  by  integrating  the  strain  between  the  point  of 
interest,  £ , and  C]>  viz. 


(I  + r~p 


(l  + G.)- 
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We  would  expect,  then,  that  the  displacements  would  approach  a 
value  given  by  the  sum 


where  the  time  variable  equals  zero  at  the  shock  arrival  time, 

. Equation  (7)  was  plotted  in  Fig.  50  corresponding  to  the 
points  of  interest.  It  should  be  noted  that  the  medium  responds 
quite  rapidly  to  this  equilibrium  condition.  The  oscillations 
will  persist  for  some  time,  which  is  as  yet  undetermined  Further- 
more, the  maximum  excursion  from  this  equilibrium  or  reference 
displacement  may  not  occur  during  the  first  four  zones  which  were 
examined . 

III.  TWO-DIMENSIONAL  QUASI-STEADY  STRESS  WAVE  PROPAGATION 

The  importance  of  the  multidimensional  character  of  the 
stress  wave  propagation  is  well  established,  however,  in  many 
cases  the  dependence  upon  some  of  the  spatial  variables  is 
relatively  weak.  Much  can  be  gained  then,  in  simplifying  the 
analysis  and  computational  effort  by  dropping  these  unimportant 
variables.  The  superseismic  region  of  the  air- induced  ground 
shock  problem  for  large  weapons  is  one  such  case  Since  the 
distances  from  ground  zero  are  generally  large,  the  divergence 
effects  (both  horizontal  and  vertical)  are  small.  Furthermore, 
the  air  blast  pressure  waveform  changes  slowly.  Under  these 
conditions  the  complex  stress  wave  problem  can  be  reduced  to  a 
two-dimensional  quasi-steady  problem;  that  is,  an  observer  moving 
along  with  the  wave  system  views  the  problem  as  being  dependent 
upon  two  variables,  the  depth  and  the  distances  behind  the  front. 
The  latter  is,  of  course,  related  to  time. 

The  following  analysis  is  divided  into  two  parts.  The 
first  part  deals  with  the  basic  interaction  of  the  shock  wave 
or  disturbance  front  with  a modulus  discontinuity,  and  develops 
a differential  equation  governing  the  strength  with  a hydro- 
dynamic  as  well  as  with  an  elastic  model.  The  second  part 
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presents  a basic  or  exact  solution  for  the  pressure  and  velocity 
fields  behind  the  shock  wave,  however,  it  is  limited  to  the  hydro- 
dynamic  model.  Thus  for  this  case,  both  the  horizontal  and  the 
vertical  displacements  are  known.  The  results  of  the  first  part 
permit  one  to  evaluate  the  differences  between  the  behavior  of 
an  elastic  and  a hydrodynamic  media  at  the  disturbance  front 
and  establish  a measure  of  the  importance  of  the  shear  capacity 
of  the  media  to  influence  the  stress  and  the  motion  in  the  media. 

A linear  and  reversible  stress- strain  law  was  selected  as 
in  the  previous  one-dimensional  nonsteady  analysis  The  wave 
diagram  for  the  quasi-steady  problem  is  illustrated  in  Fig.  51. 

We  can  now  examine  the  interaction  of  the  disturbance  front 
with  a modulus  discontinuity.  Figure  52  illustrates  this  inter- 
action for  both  the  hydrodynamic  and  the  elastic  models. 

For  the  hydrodynamic  model  the  shock  wave,  S.  , will  inter- 
act with  the  interface,  I,  across  which  the  bulk  modulus  changes 
from  M to  M+aM.  This  interaction  will  generate  a reflected 
shock  wave  which  will  propagate  back  up  into  Region  (1)  and 
a transmitted  wave  St  which  will  propagate  into  the  undisturbed 
region . 

The  solution  of  this  interaction  yields,  (after  letting 

uM-*0) 

,2  1 


which  when  integrated  results  in  the  following  stress  amplifi- 


cation relation: 


P 


1 dr 
7 c" 


Furthermore  the  vertical  (u)  and  horizontal  (v)  particle 
velocities  are 


where  the  subscript  ( )q  refers  to  the  surface  state. 

The  same  kind  of  analysis  was  performed  for  an  elastic 
media.  The  analysis  is  more  cumbersome  due  to  the  more  complex 
nature  of  an  elastic  material.  The  elastic  material  is  charac- 
terized by  two  constants.  We  will  use  Young's  modulus,  E and 
Poisson’s  ratio  , v . We  will  assume  that  Poisson's  ratio  is 
a constant,  along  with  the  density.  The  interaction  pattern 
for  this  case  is  shown  in  Fig.  52.  The  interaction  of  a dilata- 
tion of  P wave,  P^  with  an  interface,  I,  across  which  Young's 
modulus  varies  from  E to  E+AE  will  generate  a wave  system  con- 
sisting of  four  other  waves.  A dilatation  and  shear  wave  will 
be  transmitted  into  the  undisturbed  region  and  a dilatation  and 
shear  wave  will  be  reflected  back  up  from  the  interface.  The 
dilatation  and  shear  wave  velocities  are  functions  of  E and  v. 

The  solution  of  this  interaction  problems  yields  (after 
letting  AE  — *-  0 and  integrating): 


A2-  » 


where:  a = 


1 ~ 2.' 

1 - y 


- Poisson ‘ 3 ratio 


typical  result  for  the  elastic  case  is  presented  in  Fig.  53, 


US 


f 


It  should  be  noted  that  aQ  , which  is  the  normal  stress 
behind  the  disturbance  front  at  the  surface,  is  not  equal  to 
the  surface  pressure  pQ  as  it  was  in  the  hydrodynamic  case. 

This  is  due  to  the  presence  of  a shear  wave  at  the  initial 
loading  point.  The  strength  of  the  initial  dilatation  wave 
co/po  and  the  initial  shear  wave  tq/po  have  been  investigated. 

The  strength  of  the  dilatation  wave  rapidly  approaches  unity 
when  U/cQ  > 2.  Similarily  the  strength  of  the  shear  wave  becomes 
small  as  U/cQ  increases  above  a value  of  2. 

The  differences  between  the  results  of  the  hydrodynamic 
model  and  the  elastic  model  are  small  except  in  the  regions 
where  U/cq  is  small  and  where  E/Eq  approaches  (U/cG)2.  These 
are  extreme  regions  of  the  superseismic  problem  where  the  wave 
front  approaches  the  vertical  position  and  the  onset  of  sub- 
seismic  wave  propagation  occurs.  Furthermore,  a soil  may  not 
obey  either  one  of  the  fundamental  models. 

The  ground  motions  behind  the  disturbance  front  in  the 
superseismic  region  represent  very  important  input  data  to  the 
design  of  protective  construction.  It  would  appear  then,  that 
an  exact  solution  to  the  hydrodynamic  model  would  yield  reasonably 
accurate  information  on  ground  motions.  Since  the  governing 
equation  for  the  hydrodynamic  elastic  case,  one  can  hope  that 
such  an  exact  solution  can  be  found. 

Figure  51  presents  the  wave  diagram  and  defines  the  coor- 
dinate system.  This  coordinate  system  is  fixed  relative  to  the 
moving  disturbance,  that  is,  the  particle  velocity  in  the  un- 
disturbed region  is  U (in  the  horizontal  direction).  The  basic 
equations  which  govern  the  motion  behind  the  shock  front  are: 
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Momentum: 


U 


4e 

ay 


0 

0 


Conservation  of  Mass: 


<5u  , dv 


= 0 


Since  the  equations  are  linear  we  can  consider  the  special  case 
where  the  pressure  at  the  free  surface  (y=0)  is  constant  at  the 
initial  value  pQ  . The  solution  for  any  arbitrary  surface  load 
can  then  be  obtained  by  suitable  integration. 

The  solution  obtained  in  the  one-dimensional  nonsteady 
problem  suggests  that  a similar  form  exist  for  this  two-dimensional 
quasi-steady  case,  especially  since  the  former  is  a special  case 
of  the  latter.  Let: 

p = p (y) 

This  assumption  yields,  as  before 

P = P + Ay 
o — 

and 

M_  = ^ U + Q-. —z 

Mo  (A2'  !)  + (1  + O 

The  basic  solution  is  illustrated  in  Fig.  54,  where 
x 1 

5 ' b 

The  modulus  form  is  illustrated  in  Fig.  55.  This  form  can  be 
fitted  well  to  existing  site  data. 


88 


FIG  54  BASIC  SOLUTION 


Bulk  Modulus 


The  vertical  displacement;  being  representative  of  all  of 
the  results,  is  of  the  form. 


T 


,2 


) 


The  parametric  coefficient 
of  places  in  the  results, 
given  by: 


Ya^T 

TL  A 

The  horizontal 


appears  in  a variety 
displacement  r , is 


f-  = (1  + O (t-t  ! ) 


The  basic  solution  has  been  extended  to  an  exponential 
surface  load.  Typical  displacement  results  are  illustrated  in 
Fig.  56  and  57  for  a = o (the  basic  solution)  and  x = 1.  It 
should  be  noted  that  the  traj  ectories  of  the  particle  at  C =1 
are  not  influenced  as  much  as  the  surface  particle,  due  to  their 
rather  rapid  response.  Furthermore  the  direction  of  motion 
changes  rapidly  primarily  due  to  the  characteristic  bounce  of 
the  softer  surface  layer. 

Since  in  many  geologic  formations  the  properties  of  the 
material  at  some  considerable  depth  become  relatively  uniform, 
it  becomes  important  to  couple  the  above  basic  solution  to  a 
homogeneous  base  layer.  The  modulus  variation  in  the  surface 
layer  can  be  fitted  by  the  parameter,  b,  such  that  an  acceptable 
match  is  made.  Figure  58  illustrates  such  a fit  for  an  actual 
missile  silo  site.  This  particular  curve  has  a surface  modulus 
Mq  identical  to  the  corresponding  modulus  at  the  site.  However, 
the  matching  does  not  have  to  be  done  in  this  manner,  hence  we 
actually  have  two  parameters  with  which  to  fit  the  field  data 
(i.e.,  and  b)  . Furthermore,  we  can  also  tolerate  a modulus 
discontinuity  at  the  interface  between  the  base  and  surface 
layer. 
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The  existance  of  a homogeneous  base  layer  introduces  a 
definite  major  wave  pattern  which  divides  the  region  into 
descrete  zones.  These  are  illustrated  in  Fig.  59.  Zone  I 
corresponds  to  the  basic  solution  given.  The  solutions  for  the 
remaining  regions  have  yet  to  be  developed,  however,  the  pro- 
cedures have  already  been  established  in  the  one-dimensional 
nonsteady  analysis.  It  should  be  noted,  that  although  the 
surface  layer  must  be  treated  as  a hydrodynamic  material  one 
is  free  to  choose  any  material  behavior  for  the  base  layer 
which  yields  a relationship  between  the  particle  velocity  and 
the  stress  along  its  boundary  (£«£-[_).  If  an  elastic  base 
material  is  selected,  then  the  zones  in  this  layer  will  be 
subdivided  into  additional  zones,  one  of  which  is  influenced 
by  the  transmitted  shear  wave  system. 
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FIG.  59  MODIFIED  MODULUS  WAVE  PATTERN,  ELASTIC  BASE  LAYE 
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